%%author:sousanone
%%email:nj2000515@foxmail.com
%%date:2023/3/23



clc,clear all;

addpath 'D:\matlab2022\bin\cmi'

%addpath('D:\mat\bin\ptychography-master\utils')
%% FZA generation
lamda = 632.8e-6;                                        % 入射波长，单位mm米、, 632.8nm nm-um-mm-m
k0=2*pi/lamda;                                        
Nx =600; 
range1=2000;
focallength=30;                                         %焦距 
step1=range1/(Nx-1);                                     % 步距
yp1=-range1/2:step1:range1/2;
xp1=-range1/2:step1:range1/2;
[x1,y1]=meshgrid(xp1,yp1);                               % 划分网格
position1  = x1 + 1i * y1;                               % 某点的位置
lensletr = abs(position1);                               % 某点的半径
lensang  = angle(position1);
Nc=floor(Nx/2)+1;
TFZP= exp(-1i*pi*lensletr.^2/(lamda*focallength));
a=sign(real(TFZP));
lens1 =sign(1+sign(real(TFZP)));   



%% CMI reconstruction
lamda = 632.8e-6;
pix = 5.5e-3;
z1 = 20;
z2 = 30;
% % mask1=zeros(Nx,Nx);
% % mask1(Nx/2-200:Nx/2+200,Nx/2-200:Nx/2+200)=1;

mask1=circle_mask(Nx,400,Nx/2,Nx/2);


lens2=mask1.*lens1;
% % % 这啥也不是
% cmi_prior = lens2;%.*exp(1i*lens2/max(max(lens2)));

% % 相位型
% cmi_prior = double(mask1).*exp(1i*lens2/max(max(lens2)));

% % 振幅形
cmi_prior = lens2/max(max(lens2));%.*exp(1i*mask1);


obj=imread('cameraman.tif');
obj=imresize(obj,[Nx,Nx]);
obj1=double(obj).*mask1;
phase=imread('mla.jpg');
phase=double(rgb2gray(phase));
phase=imresize(phase,[Nx,Nx]);
phase=mask1.*phase;

obj2=double(obj1).*exp(1i*phase/max(max(phase)));

obj2=phase.*exp(1i*obj1/max(max(obj1)));
temp1=Propagate(obj2,'fresnel',pix,lamda,z1);
temp2=temp1.*double(cmi_prior);
raw=abs(Propagate(temp2,'fresnel',pix,lamda,z2));




iter = 40;
obj_g = ones(Nx,Nx);
support = mask1;
ew1 = obj_g.*support;
for i = 1:iter
    iw2 = Propagate(ew1,'fresnel',pix,lamda,z1);
%     figure,imshow(abs(iw2),[])
    ew2 = iw2.*double(cmi_prior);%+iw2.*(1-double(cmi_prior));
%     figure,imshow(abs(ew2),[])
    dp1 = Propagate(ew2,'fresnel',pix,lamda,z2);
%     figure,imshow(abs(dp1),[])
%     dp2 = sqrt(intensity).*dp1./(abs(dp1)+eps);
    dp2 = raw.*exp(1i.*angle(dp1));
%     
    new2 = Propagate(dp2,'fresnel',pix,lamda,-z2);
    new2 = new2.*double(cmi_prior)+new2.*(1-double(cmi_prior));
    dp1  = Propagate(new2,'fresnel',pix,lamda,z2);
    dp2 = raw.*exp(1i.*angle(dp1));
    new2 = Propagate(dp2,'fresnel',pix,lamda,-z2);
%     figure,imshow(abs(new2),[])l
%     niw2 = new2./rpm;
    niw2 = iw2 + conj(double(cmi_prior))./(max(max(abs(double(cmi_prior)).^2))).*(new2-ew2);
%     niw2 = iw2 + conj(double(ew1))./(max(max(abs(double(ew1)).^2))).*(new2-ew2);    
    new1 = Propagate(niw2,'fresnel',pix,lamda,-z1);
%     figure,imshow(abs(new1),[])
    ew1 = new1.*support+0.5*(new1-ew1).*(1-support);
    
    if mod(i,10)==0
        %figure,imshow(abs(ew1),[])
    end
end
subplot(221)
imshow(abs(ew1),[])
title('物体振幅')
subplot(222)
imshow(angle(ew1),[])
title('物体相位')
subplot(223)
imshow(abs(cmi_prior),[])
title('掩膜振幅')
subplot(224)
imshow(angle(cmi_prior),[])
title('掩膜相位')





